library(xtable)
library(dplyr)

# Load results data from Statistics Denmark 

load("dataan/data/mother_twin_res.rdata")
load("dataan/data/father_twin_res.rdata")

f_mother <- 
  read.table("dataan/data/f_stat_mothers.txt")
f_father <- 
  read.table("dataan/data/f_stat_fathers.txt")

# Create empty table for mothers 
main_tab_mom <- 
  matrix(NA, nrow = 8, ncol = 3)

# Assign column names
colnames(main_tab_mom) <-
  c("2009 Local", "2013 Local", "2014 European")

# Assign row names 
rownames(main_tab_mom) <- 
  c("Twinning in",
    "first parity",
    "No. of children (2sls)",
    "4",
    "Turnout if one child in first parity",
    "N",
    "Share with twins")

# Assign values
for(i in 1:3){
  # Effect estimate of twins 
  main_tab_mom[1, i] <- round(outtab_mom[[i + 1]][3,1], 3)
  # Assign standard errors
  main_tab_mom[2, i] <- paste("(",round(outtab_mom[[i + 1]][3,2], 3),")", sep ="")
  # Effect estimate of 2sls
  main_tab_mom[3, i] <- round(outtab_mom[[i + 1]][4,1], 3)
  # Assign standard errors of 2sls
  main_tab_mom[4, i] <- paste("(",round(outtab_mom[[i + 1]][4,2], 3),")", sep ="")
  # Assign turnout rate for parents with one child in first parity
  main_tab_mom[5, i] <- round(outtab_mom[[i + 1]][1,1], 3)
  # Assign N
  main_tab_mom[6, i] <- sum(outtab_mom[[5]][i,])
  # Assign share with twins
  main_tab_mom[7, i] <- round(outtab_mom[[5]][i,2]/sum(outtab_mom[[5]][i,]), 3)
  # Assign F-stat for weak instrument
  main_tab_mom[8, i] <- round(f_mother[i, 3])
}

# Set table up for fathers

main_tab_dad <- 
  matrix(NA, nrow = 8, ncol = 3)

colnames(main_tab_dad) <-
  c("2009 Local", "2013 Local", "2014 European")

rownames(main_tab_dad) <- 
  c("Twinning in",
    "first parity",
    "No. of children (2sls)",
    "4",
    "Turnout if one child in first parity",
    "N",
    "Share with twins")

for(i in 1:3){
  main_tab_dad[1, i] <- round(outtab_dad[[i + 1]][3,1], 3)
  main_tab_dad[2, i] <- paste("(",round(outtab_dad[[i + 1]][3,2], 3),")", sep ="")
  main_tab_dad[3, i] <- round(outtab_dad[[i + 1]][4,1], 3)
  main_tab_dad[4, i] <- paste("(",round(outtab_dad[[i + 1]][4,2], 3),")", sep ="")
  main_tab_dad[5, i] <- round(outtab_dad[[i + 1]][1,1], 3)
  main_tab_dad[6, i] <- sum(outtab_dad[[5]][i,])
  main_tab_dad[7, i] <- round(outtab_dad[[5]][i,2]/sum(outtab_dad[[5]][i,]), 3)
  main_tab_dad[8, i] <- round(f_father[i, 3])
  
}

## Print table
xtable(main_tab_mom)
xtable(main_tab_dad)